%=========================================================================%
% This program looks for a good start value for the SMM estimation
% by estimating the SMM function for 5000 potential parameter choices
% determined by a Sobol sequence over the parameter space
%=========================================================================%   

%Log on...
diary (sprintf('start_value_search_log_%s_eis075.txt',datestr(now,'yyyymmdd')));

%Code files (functions) are here

cd '/.../structural-model-code';
addpath ('...data/HRS MiCDA/Tables');
addpath ('...data/FSRDC/Tables');
addpath ('...data');
addpath ('.../structural-model-code/lininterp');

%-------------------------------------------------------------------------%
% Obtain data moments
%-------------------------------------------------------------------------%
% Treatment effect estimates for 56-64 y.o's starting from period 0 to 12
temp=csvread('te_emp_5664.csv',6,0,[6,0,18,1]);
TE_D=temp(1:end,1);
V_TE=temp(:,2).^2; %variances 
clear temp;

%LEHD control group employment rates
temp=csvread('Controllfp_levels.csv',1,0,[1,0,17,1]);
%Extract periods -5 to 12 
E_D=temp(:,1);
V_E=temp(:,2).^2; 
clear temp;

%Stack the moments (scale TE's by 40)
scaler = 10;
dtarget{1,1} = [E_D;scaler.*TE_D];

%Use I as wt matrix  
dtarget{2,1} = eye(length(dtarget{1,1}));

%-------------------------------------------------------------------------%
% Find start value using Sobol search
%
% There are five parameters describing preferences
% 1. sigma (EIS)/fixed as 1.5 HERE!
% 2. rho (AR(1) persistence of work disutility process)
% 3. sigma_v (AR(1) idiosyncratic variance of disutility process)
% 4. phi is the age trend slope
% 5. gamma is the labor disutility constant
%-------------------------------------------------------------------------%

%Set seed
rng(1,'twister');

%Generate Sobol sequence
temp = sobolset(4);
sobolseq = net(temp,5001); %inital value is 0 which will be dropped
sobolseq = sobolseq(2:end,:); %retain 5000 draws from the sequence
%Scale parameter space:
%Assume rho, phi, sigma_v in [0,1] for start value no scaling
%Assume gamma in  [-5,-1]
initparam = sobolseq;
eis = 1.5.*ones(length(sobolseq),1); %fix eis = 1.5
initparam = [eis initparam];
initparam(:,5) = -5 + 4.*sobolseq(:,end);
%Back transform rho, sigma_v 
initparam(:,2) = log(sobolseq(:,2).^-1 -1); %rho
initparam(:,3) = log(sobolseq(:,3)); %sigma_v

%Estimate SMM criteria for each potential parameter start value
SMM_criteria = zeros(length(initparam),1);
TEnonzero = zeros(length(initparam),1);

tic
parfor i = 1:length(SMM_criteria)
 [SMM_criteria(i), TEnonzero(i)] = smm_t(initparam(i,:),dtarget,scaler);
end
toc

%Output data from all runs to csv; sort by fit criterion
table = table(initparam,SMM_criteria,TEnonzero);
table_srt = sortrows(table,{'SMM_criteria'});
outdat = table2array(table_srt);
csvwrite('.../model-output/SMMcriteria_sobolseq_eis15_10wt.csv',outdat);
